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Abstract 

Attractive bonding interactions between molecules typically have inherent conservation laws 
which influence the statistical properties of such systems in terms of corresponding sum rules. We 
considered lattice water as an example and enunciated the consequences of the sum rule through a 
general computational procedure called "Molecular mean field" theory. Fluctuations about mean 
field are computed and many of the liquid properties have been deduced and compared with Monte 
Carlo simulation, molecular dynamics and experimental results. Large correlation lengths are seen 
to be a consequence of the sum rule in liquid phase. Long range Coulomb interactions are shown 
to have minor effects on our results. 
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I. INTRODUCTION 



Molecules in bulk interact with each other consistent with the underlying quantum me- 
chanical dynamics, for example, hydrogen bond, variety of sulphur bonds etc. These inter- 
actions, though quite strongly attractive compared to ambient thermal energy, still fluctuate 
only when the concerned molecules are very close to each other (of the order of Angstrom) . 
These molecular interactions are countable and categorizable. For concreteness, consider 
water molecules. Each water molecule has two hydrogen arms and two lone pair arms. In- 
teraction between water molecules is via these oppositely polarized bond arms and when a 
hydrogen arm of one molecule bonds with lone pair of a neighboring molecule, it results in 
"hydrogen bond" . So, for any number of water molecules (p) the number of hydrogen bonds 
(HB) and number of dangling bonds (DB) i.e., lone-pair and hydrogen arms which are not 
hydrogen-bonded, satisfy a "sum rule" . It states that the arms of water molecule are either 
hydrogen-bonded or not i.e., 

DB + 2HB = 4p (1) 

Now, if we consider a bulk of water molecules the above equation still holds when DB, HB 
and p are appropriately defined per unit volume. In other words, the local topology of molec- 
ular interactions implies a sum rule which is also true in the bulk for any thermodynamic 
conditions such as temperature, pressure, volume. Furthermore, this is also independent of 
other interactions in the dynamical system such as van der Waals' (vdW), Coulombic etc. 
These facts are not surprising since Eq.fll]) is a topological constraint which is insensitive to 
many details of the dynamics. 

In most descriptions of solid phase (ice or glasses), by design, the above sum rule is 
maintained explicitly. In the gas phase where the density of hydrogen bonds is negligible 
the sum rule is again trivially satisfied. If there are significant HB then the system can 
be analyzed by introducing appropriate density of dimers which still satisfy the sum rule 
thermodynamically. 



Liquid phase is typically studied using molecu 



simulations of popular models 
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ar dynamics (MD) and Monte Carlo (MC) 
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401 ] where in again the sum rule is 



maintained at every step of the dynamics and many bulk and ot 



rer correlation properties 
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Liquid phase is also described phenomenologically by a mean field which invariably is 
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taken to be density field and fluctuations there of are then postulated to obey Ornstein- 
Zernike (OZ) equation . fisj] - Furthermore, all molecular interactions have the property 
that they cannot approach each other closer than certain minimum distance, usually incor- 
porated as hard sphere repulsion in the interaction potential between the molecules. This 
in the liquid phase gives us the well known vdW term in the equation of state, often re- 
ferred to as density saturation effect [12]. Phenomenologically, such saturation effects can 
be incorporated in the OZ equation and thereby understand important components of den- 
sity fluctuations namely coordination shell peaks using approximations like Percus-Yevick, 



hypernetted-chain etc 



121 ] . But these descriptions are not as successful as MD simulations, 



primarily because by design in the construction of phenemonological models one does not 
envisage fluctuations in HB or DB densities. Consequently, they cannot satisfy the sum rule 
in large subvolumes where surface effects are negligible. The sum rule argument suggests 
that to describe water we need atleast two fluctuating fields to be self-consistent. 

Some important developments in the past on systems with hydrogen-bond' like interac- 
tions are discussed below. Wertheim's work on statistical thermodynamics of fluids with 
directional interactions consists of systematic graph-theoretic expansions of n-mer tree-like 



configurations and provides OZ- 
density fields of the system [42 



ike operator equation in terms of total density and singlet- 



44j . Another significant work on associating fluids was by 



Andersen who provided a Mayer-like density expansion technique to deal with directional in- 
teractions [l| ■ Both the descriptions work very well at high temperatures and also at densities 
corresponding to liquid- vapor coexistence. An extension to Wertheim's theory was proposed 
as SAFT 9I. I25I] wherein the equation of state is phenomenologically extended to include 
contributions from complex geometries of constituent molecules which interact through di- 
rectional forces. Another phenomenological description aimed at describing mainly the low 
temperature anamolies in fluids with characterstic orientation-dependent interactions was 
given by Truskett et al {39}] . Here in again, the equation of state is approximated with 
vdW like terms which are then phenomenologically supplemented with hydrogen-bonding 
contribution that captures certain high-density features. 

In this paper we propose to study a simple model Hamiltonian which incorporates the 
molecular hydrogen-bonding interactions. To accomodate the hard sphere repulsion we 
envisage the system on a hypercubic lattice and the model is essentially a slight generalization 
of the Pauling's model of water (3^]. The partition function corresponding to the lattice 
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model is analyzed by introducing appropriate discrete lattice fields. It is shown that the 
sum rule is automatically true in the bulk. 'Molecular mean field' (MMF) approximation 
extremizes the partition function in terms of the defined auxiliary fields. In addition, all the 
observables such as p, HB and DB densities are functionals of these auxiliary fields. One of 
the mean field equation which implies sum rule also implicates the "equation of network" 
i.e., a relation between equilibrium densities of HB and p. We study the equation of state 
and various mean field fluctuations in terms of correlation functions. We also considered long 
range Coulomb interaction and studied its consequences. Subsequently, an MC simulation 
study of the model is pursued and compared with the mean field theory quantitatively. We 
also discussed results of our analysis in the context of experiments and MD simulations. A 
brief discussion on future directions is presented at the end. 

II. MODEL FOR WATER 

On a three dimensional hypercubic lattice we define W(x) = {0, 1} corresponding to 
water being absent or present respectively, at the site x. At each occupied site we define six 
arms H a (x) = {0, ±1}, where a = {±1, ±2, ±3} corresponding to six directions around the 
site. H a (x) = corresponds to no arm on the corresponding link, +1 to that of hydrogen 
arm (H), —1 for lone pair arm (L). The constraints between W(x) and H a (x) being, 

Y,H 2 a {x) = 4W(x) (2) 

a 

H a (x) = (3) 

a 

which imply that every water has two H arms and two L arms only. A hydrogen bond is 
realized when two water molecules two lattice units apart have their H and L arms meet at 
a site, as shown in Fig.([T]). 

At any site and on any link on the lattice the possible configurations are as shown in Fig. 
(PQ). We assign an energy — v for every unpaired arm (H or L), and —A for every hydrogen 
bond. We disallow for any two H (or L) arms to meet at the same site; also, more than two 
arms are disallowed to meet at a site as in Fig.(j2J). 

To implement the constraints as shown in Fig. fl2]) in our analysis it is useful to define two 
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discrete integer fields b(x), q(x): 

b(x) = ^H 2 a (x + e a ) (4a) 

a 

= ^2H a (x + e a ) (4b) 

a 

where e a is a lattice unit vector in the direction a with the property e a = — e^. The discrete 
field b(x) counts the number of non-zero arms (approaching arms) in the neighborhood of 
site x, while q(x) measures the charge i.e, the difference between number of H arms and L 
arms meeting at site x. By construction, b(x) varies between and 6 on a three dimensional 
hyper-cubic lattice and by imposing the condition that b(x) < 2 in our analysis we ensured 
that no more than two arms can meet at a site. q(x) in turn varies between — b(x) to b(x). 
In terms of these variables the Hamiltonian is : 

H = E !) " ^( 6 (z), 2)<5(<z(a0, 0)) (5) 

X 

with constraints, 

W(x)b(x) = (6a) 
< b(x) < 2 (6b) 
q(x) = {-b(x), -b(x) +2,..., b{x) - 2, b{x)} (6c) 

The range of q(x) follows from Eq.fll]). The Kronecker delta function denoted here as S(p, q) 
is defined as 5(p, q) = 1 for p = q and otherwise. 

We now write the grand canonical partition function for our system at a finite chemical 
potential \i for water and inverse temperature (3: 



n yi ?2 \ 

x W(x),H a (x) v ' ri{x),<t>{x) 

b(x),q(x) 



— y 



exp [ KK X ), !) + & S(b(x), 2) 8(q(x), 0) + /3fiW(x) 



i—t}{x) I y^^H^x + e a ) - b(x) \ + i^z(j){x) 
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^2,H a (x + e a ) - q(x) 
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(2N + iy 
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r/(x),(f>(x) 



Y\ Z s i te (x) 



(7a) 
(7b) 
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where we have introduced two auxiliary fields r](x) and 4>(x) to impose the constraints Eq.(j4|). 
The discrete fields 77 and <fi take integer values in the range [-N, N] at every site, where N is 
any suitable large integer (greater than 8). The prime over summation refers to sum being 
restricted to constraints Eqs.(j2l [31 ED- The introduction of auxiliary fields T)(x) and (f>(x) 
allows summation over other discrete fields W(x), H a (x),b(x),q(x) within their respective 
allowed range at each site x without any restrictions from the neighborhood configurations 
i.e., as if a single site functional Z S i te . 

The Z S it e expression thus obtained is stated below. For brevity in the following expression 
jfT) is written as 77, -^0 as and further using fugacity definitions v = exp(/3z>), A = exp(/3A) 
and p = exp(/3/i), 

Z slte = l + 2ue -^cos^x) + Ae (8) 
+p {exp [ irj(x + d) + irj(x + ei) + ir}(x + e 2 ) + ir](x + e 3 )] 

(2cos(Vi0 + Vi0 - V 2 - V 3 0) + 4cos(Vi0 - Vi0)cos(V 2 - V 3 0)) 
+ . . . other orientations} 

Various terms in Eq.flS]) follow from the fact that at any site x there are only following 
contributions to the partition function : (i) unity for vacuum, (ii) v term for unpaired H or 
L arms (dangling bond), (iii) A term for hydrogen bond and (iv) fi term for water with all 
its possible orientations suitably weig hted. For the // term there are [J] arm orientations of 
having four non-zero arms around any water site and each of these in turn has [ 2 ] possible 
charge orientations corresponding to Fig.(|3]). Explicitly in the Eq.flH]) we have summed up 
these [ 2 ] charge orientations. Similarly other arm orientations are implicated in the Eq.flBJ). 
It is useful to note the following identity in implementing this computation. 

<P{x)H a {x + e a ) = J2 <P( X + e a )H a {x) = ^ 0(x) V a H a {x) (9) 

x,a x,ce x,a 

where V a /(x) = (f(x + e a ) — f(x)). The last equation follows from the constraint Eq.(j3]). 
The observables in the system such as water density p, hydrogen bond density (HB) and 
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dangling bond (DB) are calculated as follows : 





P 
HB 

DB 



A 



d 



dX \V 



1 



d 



InZ 



(10a) 
(10b) 
(10c) 



dv \V 

where V is the volume i.e., total number of lattice points. 

r]{x) and <f)(x) are discrete fields varying in the range [— N, N]. By construction the 
partition function is independent of N for N > 8. In practice it is convenient to evaluate this 
partition function by taking TV — > oo, whereupon effective r](x) and <p(x) become continuous 
fields. We implement this limiting procedure and check if the sum rule is obeyed. In the 
N — > oo limit, summation over 77 and is replaced by integrals. The resulting functional 
integral has the following trivial property : 



n 



dr](x) d(j)(x) 
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2tt 
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5^ dr](y) -[I Z 



site 







(11) 



Taking derivatives explicitly in the above gives terms proportional to u, A, /i. Since these 
are being summed at all sites each of these terms can be regrouped in terms of derivatives 
of z/, A, and /1 as : 



dv dX dfi 







(12) 



The /j, dependent term in the Eq. ffl2|) has contributions from the four neighboring sites. 
Since all sites are being summed over, the /i term in Eq.t flTT) gets each contribution four 
times. We notice that this equation is precisely the "sum rule" constraint Eq.flTJ). This 
demonstrates that the sum rule in terms of continuous auxiliary fields is automatically true. 



III. MOLECULAR MEAN FIELD THEORY 



Next, we evaluate the functional integral within the MMF approximation. The partition 
function's integrand can be envisaged as a product of field dependent phase factors at each 
site. When we enumerate them site-by-site the phase factors cancel exactly correspond- 
ing to physically allowed configurations. Evaluating along this procedure amounts to the 
standard high temperature or Mayer-like expansion. Instead we attempt an approximate 
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method wherein we first notice that if we relax the constraints Eq.(T4]) the integrand still 
peaks for the same configurations that obey Eq.(T4]) strictly. Hence in the thermodynamic 
limit, approximating the integrand suitably around the peaking configurations we may reli- 
ably estimate the partition function. This reliability can be self consistently established by 
computing the variance or correlation functions. 

The leading contribution to the functional is expected to come from the extremum which 
maximizes the integrand Z S i te - Furthermore, since we are choosing to describe fluid phase of 
the model we seek such spatial configurations which are discrete translational and rotational 
invariant in auxiliary fields. The Z site over a space-independent field configuration 77, <p is 
given by : 

Zsite\ n=n ^ = (l + 2ve -^cos0 + Ae + 90//e 4 ^) (13) 

It is evident that the maximum for Z S i te occurs at fj = = since all fugacities are 
positive. Z site at the maximum is given by Z Q : 

Z = (1 + 2u + A + 90//) (14) 

The extremization with respect to <j) is trivially true, while that with respect to fj yields, 

2v + 2A = 4(90//) (15) 

This is a consequence of the sum rule Eq.flT]) within the zeroth order approximation. Various 
densities are calculated from Eq.f llOj) with the approximate partition function Z Q . Further- 
more, using Eq. ( TT5]) : 

DB = — (16a) 

2 + 5u + 3A v ; 

9 A 

HB = — (16b) 

2 + 5u + 3A v ; 

P= 2 + 5, + 3A (17) 
Eliminating A in Eqs. (116)17p we obtain : 

HB = 2p--^ I (l-3p) (18) 

This is the "equation of network" which indeed is a different way of writing the sum rule 

Eq.([T]) in terms of dangling-bond fugacity v. Z is given by : 

5 3 \ l + i/ 1 . . 

1 + 2 y+ 2 A )-T^= l-(5p-HB) (19 > 



Rewriting it as the last term in Eq.( TT9|) shows that it is also inverse of the density of 
void sites (every water molecule in our model necessarily occupies five sites and for every 
hydrogen-bond one site is double-counted and hence is compensated in the expression). 

From the sum rule it follows < HB < 2p. Consequently, p here varies between 
and ~. The upper bound on p is indeed the highest possible density in the model while the 
lower bound is a consequence of MMF approximation, meaning that this description is self 
consistent only for densities greater than . Without loosing any generality, we choose 
A = 1 or A = e 13 i.e, measuring all energies in the units of hydrogen-bond energy. Then we 
make the observation that if temperature (4) is always positive, from Eq.([T7|) we can show 
that p is greater than ^. Furthermore, as (3 — > oo, from Eqs. (I16|ll7j) we see that p — > h, 
HB — > | and DB — > 0, showing that every water molecule's arms are all hydrogen-bonded 
leaving no void sites. This saturation at p — | is verified to be exactly true by explicit 
construction of such configurations. We find that the highest density is not that of a unique 
crystal configuration. Instead, there are infinitely many configurations corresponding to 
different spatial and orientational arrangements of water molecules. 

The equation of state in terms of densities is given by : 

Vo ~~P HZo} ~ ln(l - 5p + HB) - ln(HB) (20) 

where V is the pressure and (3 is written in terms of natural energy units of the model. 

By considering the case where H and L arms are neither energetically favored or sup- 
pressed i.e., v = and v = 1, the entropy per site is given by, 

In the limit (3 — > oo, p reaches the maximum value | and the entropy at the highest density 
tends to a constant value ln(|). The constant | compares exactly with that of Pauling's 



estimate for tetrahedral ice i.e. L 2 2 ^ = 1.5 

1 1 lo 



30| and agrees well with the numerical estimate 



by Nagle i.e., 1.50685±0. 00015 [26]. We note that these results are independent of dimension, 
hence in two dimensions it also compares with the exact result of Lieb on 'square ice' 21 1. 
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A. Fluctuations 



Next, we evaluate the functional integral of the partition function by considering small 
fluctuations around the mean field and approximating them by a Gaussian i.e., one-loop 
correction Q. Expanding Z site up to quadratic terms, we obtain : 



■'site 



I : 2r + A 90// \ 2i>[ -w?-y-y) + X (-2tT] - 2r] 2 ) 



+90/i 



,,77 -877 



4z77 + ~i ^ - ( Yl ^ V 

a \ a 

4fev^) 2 -|E(v^ + ^(Ev t . 

\ a / a \ a 

Z a exp i -Z77(x) (V + 2A' - 4/x' j - ^y- ^ V a ?7 

(V + 4A' + 16//) r,\x) + ^- ( V„r, 



a 



15 



16/i 



2" 



(22) 



a \ at 

where v = 2u/Z a , X' = X/Z a , \i = 90fi/Z o are the reduced fugacities; such that all of them 
are less than 1. 

Inserting the above expression for Z site in Eq.(j7j) and evaluating the resulting Gaussian 
integral by Fourier transform in a periodic box the pressure V is given by : 



V = \\ 1 < Z °) - \ j M^(A)) + ln(P^(A))] 



where 



P vri( A ) 



2/9 A „ , ' . / 9 / ^ A 



96/2 



(23) 



(24) 



+1/ + 4A +16/2 - ( 1/ + 2A - 4/i 



A (1 - A) + v 



(25) 
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and A = - > (1 — cos(fci)). Thereupon the densities p, HB, DB are calculated using Eq. ffTO]) . 
6 z — ' 

i=l 

We compare these results with that of MC simulations in a later section (see sec JIV Al on 
MC). 



B. Correlation functions 



The position space correlation functions for rj(x) and <p(x) fluctuations i.e., G m {r) and 
G^r), to the leading order, are given by : 



G m (r) = <77(0)77(r)> 
G H {r) = (0(O)0(r)> 



Ak-r 



(27T)3p (fc) 



(2vr) 3 p^(fc) 



(26) 
(27) 



We note that to zeroth order // ~ p, A' ~ HB, z/ ~ DB and hence, using Eq.( fl5|) : 



P w (fc) ~ 64p 



96p 
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A(l-A) 



A 



20(^-p)7 ' 8(^-p)2 



+ 



5(DB) 
96p 



(28) 
(29) 



The asymptotic behavior for large r of G m and G^ correlators can be obtained by 
pursuing small- expansion of the integrand and noting that for small k, A ~ i ^\ A; 2 . The 
G w correlator for large r is, 



G m (r) oc sin(o;ir) 



where 



mi 



sin I -tan 1 W — I — — p 



'50 / 9 



10 " 



27 V25 



15 -P 



cos -tan 



'50 / 9 



27 V25 



P 



(30) 

(31) 
(32) 



This shows G m has periodic peaks reminiscent of the coordination shell peaks whose ampli- 
tudes are exponentially falling off. 
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The correlator takes the following asymptotic form for large r, in addition to oscil- 
latory behavior prominent at short distances : 



Gssir) oc ; m 2 



\ 




(33) 



We make the observation from Eqs. ( l28|29i) that the correlation functions Eqs. ( 1261271) of the 
mean field theory diverge if p — > i.e., even the local fluctuations about the mean field 
are very large rendering the approximation invalid. Indeed the theory fails well before that 
because it violates the sum-rule already at p — • This is the reason why our mean field 
theory fails for the gas phase and thus does not describe any liquid-gas transition. To sum 
it, in our model MMF approximation shows density saturation and the description is more 
accurate at higher densities only. 

Next, we compute some important physical correlations such as water- water and charge- 
charge. 

The non-zero value of VT-field at each site x, y picks only the term proportional to the 
fugacity p in Z site expression Eq.flB]) at both sites. The {. . .} include same set of terms as in 
Eq.flH]). We now make the one- loop approximation of retaining terms only upto quadratic 
in fields and further, subtract out single-site averages i.e., (W(x)W(y)) — {W (x)) {W (y)) . 
Thereby, correlation of water fluctuations to leading order is given by : 

— 7T 

A similar computation for the q(x) field can be done, wherein the non-zero q value picks 
terms proportional to the fugacity v in Z site at either sites (EqJHJ). To the leading order, 

7T ^ 

/^3l p ik-f 
W? TJf) (36) 

— 7T 

Note that the single-site average over q(x) field is zero. We remark that in this model W{x) 
being charge neutral gets contribution from the neutral r](x) field, while q(x) from that of 
the charged <f>(x) field. 
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C. Coulomb interaction 

In this section we consider the influence of long range Coulomb interaction potential 
between the charges H, L given in terms of charge e as : 



where prime over summation means x 7^ y. In our model we envisage the charges at the tip 
of water's H or L arms. This can be incorporated in our analysis by using an auxiliary field 
technique : 



where the laplacian operator nx(x) = J2 a (x( x + e ") ~ x( x )) — J2 a ^aX( x ) an d m is a 
parameter that regulates the range of interaction. The interaction potential behaves as ^—1 
for large distances r, which also reduces to Coulomb interaction when m = 0. If the lattice 
constant and m are both taken to be zero, then it reduces to exact Coulomb interaction for 
all r. 

By inserting the above in our partition function Eq. ([7]) all the interactions of water degrees 
of freedom remain unchanged with the following transformation 77 — > 77, (fi — > (fi + \\/ fie 2 - 
The extremum of the new partition function is still at fj = <fi = \ = 0. The leading zeroth 
order term remains unchanged; the one-loop correction about the mean field gets additional 
contributions due to quadratic terms corresponding to (fix an d XX i n the Gaussian expansion, 
given by : 




(37) 




(38) 




x _ a 




(39) 



(40) 



The pressure V with Coulomb interactions is, 



V 
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where 



P(A) 



p p p2 

r <t><t> r XX r <f>x 



96/i' 



12A + m 2 



A(l-A) + 



5z/ 



/3eV 



96/i' 12A + m 2 



A(l-A) 



which to the leading order can be written as : 



A/1 AN 5(DB) / /3e 2 (DB) . x , x 



p(£) ~ ^(i2A + m 2 ) _.. .... . 

y ' 5 V L 96p V 12A + ^ 2 

The two-point correlations of the fields in momentum space are given by : 



^$i)x$ 2 )) 



P{h) 

P^jkJSjh + k 2 ) 

P(h) 
P 0x (fc 1 ) ( 5(fc 1 + k 2 ) 

P(h) 



(42) 



(43) 



(44a) 
(44b) 
(44c) 



The factor e 2 being small the strength of Coulomb interactions in comparison with other 
interactions is weak. Furthermore, all charge effects are proportional to /3e 2 (DB) wherein 
/3(DB) is always finite since as /3 — > oo, DB — > 0. Consequently, effect of Coulomb in- 
teractions on thermodynamic properties such as equation of network and equation of state 
is small. Correlation functions are not perturbed very much at large distances. It does, 
however, shift the coordination peaks slightly. 



IV. MONTE CARLO SIMULATION 

The hydrogen-bond network model with all the accompaning constraints is simulated in 
the background of a cubic lattice using standard Monte Carlo (MC) procedure for grand 
canonical ensemble i.e. jlVT ensemble 10|, |20). The moves implemented at each simulation 
step on a randomly chosen site are : 

• if there is no water, "insert" one provided there are four arms free around the site 

• if there is water, with equal probability either "delete" water or "rotate" the arms to 
another possible allowed configuration as dictated by Eq.([5]). 



14 



The new configuration is accepted with a probabilit y of Boltzmann weight over energy change 



in compliance with detailed balance condition |10l . |45|. All energies /3, z>, \x are calculated in 
units of hydrogen-bond energy i.e., A = 1. Furthermore, we considered specifically v = case 
thereby making the theory essentially parameter-free other than temperature and chemical 
potential which are varied in steps as per needs of the simulation. 

The thermodynamic observables of the fiVT ensemble i.e., number of water molecules 
p, number of hydrogen bonds HB, total energy E are obtained as averages from the sim- 
ulation. A significant equilibriation time is allowed before the production run begins. In 
the production run the configurations are sampled about every 50 — 300 MC steps over a 
simulation time of 10 5 — 10 6 MC steps in order to carry out averaging procedure. Running 
averages computed at every sampling step allow determination of the efficiency of sampling 
procedure over the extent of simulation time. It is ensured that there is no observable overall 
rise or fall in the averages and variances and that the results smoothly converge to within a 
relative error of 10 — 2 — 10 3 . Then the same are noted down for record. At every step various 
configuration' statistics are recorded such as number of molecules with one, two, three, and 
four hydrogen bonds. The on-the-run variances in all the quantities are computed using 
standard discrete sum formula. 

The size dependence of the averages is ascertained and an optimal lattice size of 20 sites 
per side is found to be closely reproducing averages upto fourth decimal relative to bigger 
lattice sizes. 

A 'successive seeding' procedure is employed to adiabatically progress in the parameter 
space of the model. For example, to carry out simulations at a fixed temperature and a 
range of chemical potentials an initial run carried out at low densities and equilibriated over 
long time is used as a 'seed' i.e., input configuration for the successive run. This is repeated 
to reach higher densities. This procedure accelerates equilibration considerably compared 
to any random seed configuration. Furthermore, since a range of chemical potentials are ex- 
plored the step size in the parameter value is appropriately adjusted such that a quench-like 
situation is avoided. Thereby, the system evolves smoothly in configuration space without 
any unwanted domains persisting. The end averages remain unchanged when the seeding 
procedure is carried out in an alternative parameter space, for example, instead of chemical 
potential temperature can be varied in small steps. This confirms the absence of any pos- 
sible bias created by our successive seeding procedure in most part of the parameter space 
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(except near first order phase transitions where hysteresis exists). 

In order to obtain equation of state using MC we employed fixed temperature simulations, 



fol 
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owea 



28 



by application of Gibbs-Duhem procedure to obtain pressure versus density curves 



291 ] . In the pVT ensemble temperature is fixed and a range of chemical potentials 
are explored starting from zero density where pressure can be normalized to the value zero. 
Employing the successive seeding procedure we obtain various thermodynamic averages as 
functions of chemical potential. The Gibbs-Duhem procedure involves integrating the p vs. 
p curve so that pressure can be obtained from the following relation : 

Pi 

V = f p(p) dp (45) 



where the pi corresponds to p = and V = and jlf to that of the desired p. A set of 
representative p vs. p curves are shown in Fig. (J3J). Each such curve has two prominent 
shoulders where considerable slope change occurs. At very high p values no more equilib- 
riated configurations could be traced and the density starts shooting up to the saturation 
value. 

We also made preliminary investigation in studying phase transitions in the model. As 
shown in Fig. (01) the model exhibits discontinuity in p(p) for temperatures /3 > 2 (in units 
of hydrogen-bond energy). For instance, at /3 — 3 the system jumps from a low dense 
(p ~ 0.025, h ~ 1.2) to higher density (p ~ 0.16, h ~ 1.7) where h = — is average number 
of hydrogen-bonds per molecule. Hysteresis is also observed when the system progresses 
first from low-dense to high-dense state upon increasing p and then retracing the path by 
decreasing p. This indicates the presence of first-order phase transition in the region. We 
interpret this as liquid-gas transition in the model. 

A. MMF and MC 

One of the important expositions of the MMF theory is the "equation of network", a 
functional relationship between total water density and HB density. From Eq. ffT8~]) . at v = 1, 
the mean field equation of network is given by, 

HB = ^ - \ for p > 1 (46) 
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In MC simulation we identify the liquid-gas transition for temperatures /3 > 2. As shown 
in Fig.(jH), for = 3.0 at p, ~ —1.91 water density jumps from 0.025 to 0.16. The jump gap 
increases with /?, so does the higher density state to which system jumps. Within our limited 
exploration of the phase diagram we identify that water densities p > 0.16 correspond to 
liquid phase, only above which MMF theory is self-consistent (Eq j4"6"j) . Furthermore, the 
linear relationship between p and HB is accurate at high pressures as confirmed by MC 
simulation. As noted earlier, MMF approximation complies with the sum rule of the system; 
the same is exactly satisfied by MC simulation at every configurational move and hence over 
ensemble average. Thus, the agreement between mean field (zeroth order and one loop 
correction) and MC over the equation of network is inescapable. However, since the mean 
field description holds forte in thermodynamic limit, the HB density in the above equation 
corresponds to lowest free energy (or high pressure) states for each water density. 

Next, we compare the MMF computation of equation of state with that of MC simulation. 
From theory, at v = 1, 



The comparison is put forth in Fig.([7|). It shows that the locus of high pressure states at 
each density in MC simulation has the same profile as in MMF theory (both in zeroth order 
and one- loop'). Further, the qualitative agreement is good only in the high density regime 
where the equation of network too is accurate. A quantitative comparison of equation of 
state between MMF and MC is unreliable because MMF pressure absurdly vanishes at p — | 
whereas, physically, the pressure is zero in this model only at p = 0. As discussed earlier, 
MMF approximation fails for small densities. Therefore, a consistent normalization between 
various schemes of calculation is not present. Thus, the qualitative picture obtained from 
MMF calculation is only indicative, nevertheless consistent with MC results. 

The spatial correlation functions are also computed from the MC simulation and com- 
pared with that of analytical expressions obtained within the mean field approximation. 
Firstly, we note that the underlying hyper-cubic lattice dominates all correlations especially 
at short distances. Then, to extract the rotational invariant part of any function f(x) we 



o 




for p > - 



(47) 
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implement the following projection procedure : 

R(r) = © (\x\ - r) 6 ((r + 5) - \x\) (48a) 

X 

= Wr)Y, © (N " r) 6 ((r + 5) - |x|) (48b) 

where |x| = \/$^ xf and 9 is Heaviside step function defined as 9 (x — a) = 1 for x > a and 
for x < a. The important correlations in the model are < W(0)W(r) > and < q(0)q(r) >. 
The /c-space integrals in Eqs. fl35| [361) are computed numerically and the MMF correlation 
functions are compared with those of MC simulation. The positions of coordination peaks 
in both are in agreement. In addition, charge correlations show similar exponential fall-offs 
asymptotically. Fluctuations are seen to increase as we approach saturation density in both 
cases, however they are higher in simulation. 



V. MMF, EXPERIMENTS AND MD 



MMF theory predicts correlation length of charge fluctuations i.e., in (qq) correlation 
function Eqs.( l33|l3T)j) . These fluctuations could give rise to long distance correlations in 



other charged fields in the system such as in dipole, as seen in MD simulations (wherein two 
correlation lengths of order 5.2 A and 24 A have been observed, former being 10 times larger 
than the latter in strength) 17|. To relate to MMF we make the following observations. 
Water molecule on average participates in 3.6 hydrogen-bonds in ambient conditions 14, 
24) i.e., — = 0.4. Using this value in 1712 of Eq. (1331) the correlation length turns out 
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to be approximately 2.02 x lattice constant in MMF. We next estimate the lattice constant 
in physical units by comparing the position of highest coordination peak in (W(0)W(r)) 
correlation obtained from MMF to that of experimental data. This procedure gives a lattice 
constant equal to -p^r A. Hence in physical units MMF predicts a correlation length of 
about 3.2 A for liquid water. It should be noted that these predictions are not robust as the 
coefficients such as ^ in Eq. (l33p might vary with topology of the underlying lattice. 

We also pursued a preliminary comparison of the equation of network from experiments 
and MD simulations. The density of hydrogen bonds is indirectly probed and inferred in 



various experiments 



aa 
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varying external conditions 
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and also computed in MD and MC simulations under 
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24j. The latter works used different definitions for 
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hydrogen-bond computation such as energy-based, geometry-based or hybrid. The data is 
put in the perspective by converting the mass densities to number densities using known 
radius of a water molecule. As shown in Fig.(jS]), in the region of high molecular density 
i.e., corresponding to liquid water we find that HB density and water density are linearly 
related to each other. A linear fit function is used for HB vs. p curve and compared with 
that of the MMF equation Eq. (1181) . We find that the dangling bond fugacity v is in the 
range (0.06,0.18), implying that the corresponding energy v is positive and large compared 
to thermal energy i.e., dangling bonds are highly disfavored in liquid water. 

Another useful quantity namely fraction of water molecules with i hydrogen bonds can 
also be calculated. Consider a water molecule in a configuration where its z-number of arms 
are hydrogen-bonded to neighboring molecules and its other (4 — i) arms remain dangling 
type. A weight can be associated with each such configuration defined in terms of appropriate 
site fields and summed over all possible orientations of the molecule. We denote this weight 
averaged with respect to the full partition function for each i as pi. For instance, the averaged 
weight assigned to a molecule which is hydrogen-bonded to only two other molecules is given 
by : 

P2 = Yl 6 ( b ( x + e °i)' °) 6 ( b ( x + e ^)> °) 6 ( b ( x + e «s)> 2 ) + e <* 3 )> °) 

{a>i,a2,..;Ot6} 

5(b(x + e a4 ), 2) 8(q(x + e a J, 0) 5(b(x + e as ), 1) 8(b(x + e ae ), 1)) (49) 

The prime over summation means dissimilar a. The probability for a z-bonded molecule 
at any site x is probability that any two directions around central site have zero arms, 
each denoted by 5(b(x + e a ), 0), that i other directions have a hydrogen-bond denoted by 
S(b(x + e a ), 2) 8(q(x + e a ), 0) and that remaining (4 — i) sites are of dangling bond denoted 
by S(b(x + e a ), 1). The summation over the set {ati, a 2 , . . . , a 6 } implies summing over all 
possible rearrangements of hydrogen-bonds, dangling bonds among all the directions. With 
[2] ways of choosing empty site, LI ways for two hydrogen bond sites the P2, to the leading 
order, is given by : 

(HB) 2 (DB) 2 = 15(6)p 4 /i 2 (2 - h) 2 (50) 
where h = — is average number of hydrogen bonds per molecule. Similarly, other p^s can 
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be enumerated and computed upto leading order. For i = 0, . . . , 4, 

Pi ~ 15 4 p 4 ^(2 - hf- 1 



(51) 



Thereupon, /j which is fraction of i-bonded molecules can be calculated from the relation, 



fi 



V Pi 



h l (2-h) 



4-i 



(52) 



Note that the above expression is obtained to zeroth order approximation within the model. 
There exist one-loop corrections to it that can be calculated from the MMF theory, but 
numerically they are small. These distributions agree very well with MD simulations j^J. 
Further, probabilities for various cluster configurations (n-mers) can also be calculated within 
the MMF theory along the same lines as above. Liquid water is known to form molecular 
clusters such as trimers, tetramers, pentamers [33j. 



VI. DISCUSSION 



We analyzed the statistical mechanics of molecules with attractive bonding interactions. 
We argued that these systems typically have sum rules built-in which have to be respected 
in the analysis. In addition, there are other constraints, such as those shown in Fig.(T5]) for 
water, all of which are captured by defining fields such as b(x) and q{x) as in Eq.(TJ|. The 
analysis yields ultimately a functional integral in terms of auxiliary fields r)(x), <f>(x) which 
are conjugates of the discrete fields b(x) and q(x) respectively. A simple-minded mean field 
analysis of the auxiliary fields is shown to reproduce all major thermodynamic properties of 
the system, which are verified to be reasonably accurate in liquid phase and near saturation 
density. Principally, "equation of network" which is another manifestation of the "sum 
rule" plays an important role. We showed that all these consequences are known to be 
qualitatively true in experiments and MD simulations as well. 

The MMF approximation is shown to be the leading order in a formal Feynman loop ex- 
pansion of the functional integral with necessarily small parameters called reduced fugacities. 
Correlation functions and their properties have been computed and shown to qualitatively 
agree with MC simulation. It is however observed that many of the details in correlation 
functions do depend upon the underlying lattice, but the analytic structure is amenable to 
interpretation in the continuum as well. Pressure and correlation lengths have large con- 
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tribution from the fluctuations of DB and p in the form — . We note that any long range 
interactions such as Coulombic have very little quantitative effect on our results. 

The MMF theory is seen to be inconsistent at low densities, consequently fails to describe 
the gas phase of the model and the corresponding liquid-gas phase transition. In the explicit 
model considered here we do not have a crystalline solid phase, but we may have a glassy 
phase. Namely, at zero temperature the lowest energy configuration is infinitely degenerate. 
In this region MMF theory does not show any liquid-glass first order phase transition, 
however it does indicate a second order transition as seen from q(x) correlator. In Eq. fj33jl 
for (3 — » oo: p — > |, DB — > implying m 2 — > 0. In MC simulation we did observe dynamical 
slowing down as we reach saturation density, perhaps indicative of a phase transition. 

The MMF technique can be applied to associating liquids which may constitute variety 
of molecules and hence with variety of sum rules. Therefore, we expect a variety of dangling 
bond fluctuations contributing to pressure and correlation lengths. Indeed, there may even 
be variety of phase transitions. Essentially, the MMF technique is amenable to analysis of 
all such systems. 
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IX. FIGURES 
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FIG. 1: Allowed configurations : (grey circle) water site, (+) H arms, (-) L arms and (0) no arms; 
two types of dimer hydrogen-bonded via (+) H arm and (-) L arm; (right bottom corner) unit 
vectors' definition on cubic lattice. 











X- 






—> 












-X 




X- 


X- 





























FIG. 2: Disallowed configurations : (+) H arms (or (-) L arms) of two neighboring molecules 
meeting at a site, more than two non-zero arms meeting at a site. 
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FIG. 3: A set of orientations consistent with Eqs, (|2|3p and corresponding to representative weight 
in fx term of Eq.flSJ). 
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FIG. 4: First order phase transition in MC simulation seen for /3 > 2: isotherms are (3 = 3.0, 2.0, 1.5; 
(red) (3 = 3.0 isotherm discontinuity at fl = —1.91, density changing from p ~ 0.025 (h ~ 1.2) to 
p~ 0.16 (h ~ 1.7). 
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FIG. 5: Equation of network for v = : (green, dashed lines) MC simulation isotherms (temper- 
atures increasing from bottom to top); (blue, dotted line) mean field equation; (magenta, open 
triangles) with one-loop correction. 
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FIG. 6: Equation of network : (blue, filled circles) Experiments and (magenta, filled triangles) MD 
simulations. 
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FIG. 7: Equation of state : (green, dotted line) MMF zeroth order; (magenta, filled triangles) 
with one-loop correction; (red lines) isotherms (temperature increasing from top to bottom) in MC 
simulation; (blue, dashed line) locus of high pressure states. 
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FIG. 8: MMF : water fluctuations' correlation function 
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FIG. 9: MC : water fluctuations' correlation function 
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FIG. 10: MMF : charge fluctuations' correlation function 
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FIG. 11: MC : charge fluctuations' correlation function 
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